Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/cpp_libs/lin_alg.lha/lin_alg/LinAlg.h

C/C++ Source or Header | 1993-08-09 | 14.3 KB | 478 lines

// This may look like C code, but it is really -*- C++ -*-
/*
 ************************************************************************
 *
 *              Linear Algebra Package
 *
 * The present package implements all the basic algorithms dealing
 * with vectors, matrices, matrix columns, etc.
 * Matrix is a basic object in the package; vectors, symmetric matrices,
 * etc. being considered as a matrix of a special type.
 *
 * Elements of the matrix are arranged in memory in the COLUMN-wise
 * fashion (in the FORTRAN spirit). In fact, it makes it very easy to
 * feed the matrices to FORTRAN procedures, which implement more
 * elaborate algorithm.
 *
 * Unless otherwise specified, matrix and vector indices always start
 * with 1, spanning up to the specified limit.
 *
 ************************************************************************
 */

#pragma once
#pragma interface

#include "myenv.h"
#include <std.h>

#define REAL float            // Scalar field of the Linear Vector
                    // space

typedef REAL (*USER_DEFINED_REAL_F)(REAL);
typedef double (*USER_DEFINED_DOUBLE_F)(double);

class Vector;                // Vector over the real domain

class MatrixRow;            // A row of the matrix
class MatrixColumn;            // A column of the matrix
class MatrixDiag;            // The diagonal of the matrix
class MatrixPivoting;            // For the determinant evaluation


/*
 *------------------------------------------------------------------------
 *            Basic class of matrix
 */

class Matrix            // General type matrix
{
  friend class MatrixRow;
  friend class MatrixColumn;
  friend class MatrixDiag;
  friend class MatrixPivoting;

private:            // Private part
  int valid_code;            // Validation code
  const int MATRIX_val_code = 575767;    // Matrix validation code

protected:            // May be used in derived classes
  short nrows;                // No. of rows
  short ncols;                // No. of columns
  short row_lwb;            // Lower bound of the row index
  short col_lwb;            // Lower bound of the col index
  char * name;                // Name for the matrix
  int nelems;                // Total no of elements, nrows*ncols
  REAL ** index;            // index[i] = &values(0,i) (col index)
  REAL * elements;            // Elements themselves

  void allocate(const short nrows, const short ncols,
        const short row_lwb = 1, const short col_lwb = 1);

public:            // Public interface

                // Constructors and destructors
                    // Make a blank matrix
  Matrix(const short nrows, const short ncols);        // Index start with 1
  Matrix(const short row_lwb, const short row_upb,    // Or with low:upper
     const short col_lwb, const short col_upb);    // boundary specif.

  Matrix(const Matrix&  another);    // Make a new blank matrix a la
                    // the other one
  Matrix(const char * file_name);    // Read the matrix from the file

  ~Matrix();                // Destructor

  void set_name(const char * name);    // Set a new matrix name

                    // Erase the old matrix and create a
                    // new one according to new boundaries
  void resize_to(const short nrows, const short ncols);    // Indexation @ 1
  void resize_to(const short row_lwb, const short row_upb, // Or with low:upper
         const short col_lwb, const short col_upb);// boundary specif.
  void resize_to(const Matrix& m);            // Like other matrix m


  void is_valid() const
  { assure(valid_code == MATRIX_val_code,"Invalid matrix"); }

                // Status inquires
  int q_row_lwb() const            { return row_lwb; }
  int q_row_upb() const            { return nrows+row_lwb-1; }
  int q_nrows()    const            { return nrows; }
  int q_col_lwb() const            { return col_lwb; }
  int q_col_upb() const            { return ncols+col_lwb-1; }
  int q_ncols()    const            { return ncols; }

  int q_no_elems() const        { return nelems; }

  const char * q_name() const        { return name; }

                // Individual element manipulations
  REAL& operator () (const int rown, const int coln) const;


            // Element-wise matrix operations

                // Matrix-scalar arithmetics
                // Modify every element of the
                // Matrix according to the operation
  Matrix& operator =   (const double val);    // Assignment to all the elems
  Matrix& operator -=  (const double val);    // Add to elements
  Matrix& operator +=  (const double val);    // Take from elements
  Matrix& operator *=  (const double val);    // Multiply elements by a val

                // Comparisons
  int      operator ==  (const double val) const;
  int      operator <  (const double val) const;
  int      operator >  (const double val) const;

                // Other element-wise matrix operations
  Matrix& clear(void);            // Clear the matrix (fill out with 0)
  Matrix& abs(void);            // Take an absolute value of a matrix
  Matrix& sqr(void);            // Square each element
  Matrix& sqrt(void);            // Take a square root

  Matrix& user_function(USER_DEFINED_REAL_F uf); // Apply a user function
  Matrix& user_function(USER_DEFINED_DOUBLE_F uf);    // to each element

                // Element-wise operations on two matrices
  Matrix& operator = (const Matrix& source);    // Assignment
                    // Arithmetics
  friend Matrix& operator += (Matrix& target, const Matrix& source);
  friend Matrix& operator -= (Matrix& target, const Matrix& source);
  friend Matrix& add(Matrix& target, const double scalar,const Matrix& source);
  friend Matrix& element_mult(Matrix& target, const Matrix& source);
  friend Matrix& element_div(Matrix& target, const Matrix& source);

                    // Comparisons
  friend int  operator == (const Matrix& im1, const Matrix& im2);
  friend void compare(const Matrix& im1, const Matrix& im2,
              const char * title);
  friend void are_compatible(const Matrix& im1, const Matrix& im2);


                // True matrix operations
                // (on a matrix as a whole)

                    // Transposition
  Matrix transpose(void) const;

  Matrix& operator *= (const Matrix& source);    // Inplace multiplication
  friend Matrix operator * (const Matrix& A, const Matrix& B);

//  Matrix& operator *= (const DiagMatrix& diag);    // Multiply by diagonal matrix



                // Compute matrix norms
  double row_norm(void) const;        // MAX{ SUM{ |M(i,j)|, over j}, over i}
  double norm_inf(void) const        // Alias, shows the norm is induced
         { return row_norm(); }        //     by the vector inf-norm
  double col_norm(void) const;        // MAX{ SUM{ |M(i,j)|, over i}, over j}
  double norm_1(void) const        // Alias, shows the norm is induced
         { return col_norm(); }        //     by the vector 1-norm
  double e2_norm(void) const;        // SUM{ m(i,j)^2 }, Note it's square
                    // of the Frobenius rather than 2. norm

  friend double e2_norm(const Matrix& m1, const Matrix& m2);
                    // e2_norm(m1-m2)

  double determinant(void) const;    // Matrix must be square one

                // Make a matrices of special kind
  Matrix& unit_matrix(void);        // Matrix needn't be a square
  Matrix& hilbert_matrix(void);        // Hilb[i,j] = 1/(i+j-1); i,j=1..max

                // Row, Column operations
  MatrixRow    a_row(const int rown) const;    // Get a matrix row
  MatrixColumn    a_column(const int coln) const;    // Get a matrix column
  MatrixDiag    diag(void) const;        // Get a matrix diagonal

                // I/O: write, read, print
                      // Write to a file
                    // "| command name" is OK as a file
                    // name
  void write(const char * file_name,const char * title = "") const;
  void info(void) const;        // Print the info about the Matrix
  void print(const char * title) const;    // Print the Matrix as a table

};

/*
 *------------------------------------------------------------------------
 *            Inline Matrix operations
 */

inline Matrix::Matrix(const short no_rows, const short no_cols)
{
  allocate(no_rows,no_cols);
}

inline Matrix::Matrix(const short row_lwb, const short row_upb,
              const short col_lwb, const short col_upb)
{
  allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb);
}

                    // Make a new Matrix like the
inline Matrix::Matrix(const Matrix& old)    // old one
{
  old.is_valid();
  allocate(old.nrows,old.ncols,old.row_lwb,old.col_lwb);
}

                // Resize the matrix to accomodate to a pattern
inline void Matrix::resize_to(const Matrix& m)
{
  resize_to(m.q_row_lwb(),m.q_row_upb(),m.q_col_lwb(),m.q_col_upb());
}

inline REAL& Matrix::operator () (const int rown, const int coln) const
{
  is_valid();
  register int arown = rown-row_lwb;        // Effective indices
  register int acoln = coln-col_lwb;

  if( arown >= nrows || arown < 0 )
    _error("Row index %d is out of Matrix boundaries [%d,%d]",
       rown,row_lwb,nrows+row_lwb-1);
  if( acoln >= ncols || acoln < 0 )
    _error("Col index %d is out of Matrix boundaries [%d,%d]",
       coln,col_lwb,ncols+col_lwb-1);

  return (index[acoln])[arown];
}

inline Matrix& Matrix::clear(void)    // Clean the Matrix
{
  is_valid();
  memset(elements,0,nelems*sizeof(REAL));
  return *this;
}

inline void are_compatible(const Matrix& im1, const Matrix& im2)
{
  im1.is_valid();
  im2.is_valid();

  if( im1.nrows != im2.nrows || im1.ncols != im2.ncols ||
      im1.row_lwb != im2.row_lwb || im1.col_lwb != im2.col_lwb )
    im1.info(), im2.info(), _error("The matrices above are incompatible");
}

                // Assignment
inline Matrix& Matrix::operator = (const Matrix& source)
{
  are_compatible(*this,source);
  memcpy(elements,source.elements,nelems*sizeof(REAL));
  return *this;
}

/*
 *------------------------------------------------------------------------
 *        Friend classes - MatrixRow, MatrixCol, MatrixDiag
 */

class MatrixColumn        // Special representation of a Col of the
{                // matrix
  friend class Matrix;
  friend class Vector;

  const Matrix * matrix;        // The matrix I am row of
  REAL * ptr;                // Pointer to the a[0,i]

                // Private constructor
                  // Note there are no public constructors;
                // MatrixColumns  are always
                  // created via the Matrix::a_col()
                // function
  MatrixColumn(const Matrix* mp, const int col);

  public:
                // Note how to construct MatrixColumn of
                // the matrix 'm':
                //        m.a_column(10)

  ~MatrixColumn()    {}

                    // Assign a value to all the elements
                    // of the Matrix Col
  friend void operator  = (const MatrixColumn& mc, const int val);
                      // Modify the elements in the col
  friend void operator +=  (const MatrixColumn& mc, const int val);

  void operator = (const Vector & vec);    // Assign a vector to a matrix col
};

                // Construct the MatrixColumn
inline MatrixColumn::MatrixColumn(const Matrix * mp, const int col)
{
  matrix = mp;
  matrix->is_valid();

  register int colind = col - matrix->col_lwb;

  if( colind >= matrix->ncols || colind < 0 )
    _error("Column #%d is not within the matrix",col,((*matrix).info(),0));

  MatrixColumn::ptr = &(matrix->index[colind][0]);
}

inline MatrixColumn Matrix::a_column(const int coln) const
return c(this, coln)
{}


class MatrixRow            // Special representation of a Row of the
{                // matrix
  friend class Matrix;
  friend class Vector;

  const Matrix * matrix;        // The matrix I am row of
  const int inc;            // if ptr=@a[row,i], then
                    //    ptr+inc = @a[row,i+1]
                    // Since elements of a[] are stored
                    // col after col, inc = nrows
  REAL * ptr;                // Pointer to the a[row,0]

                // Private constructor
                  // Note there are no public constructors;
                // MatrixRows  are always
                  // created via the Matrix::a_row()
                // function
  MatrixRow(const Matrix* mp, const int row);

  public:
                // Note how to construct MatrixRow of
                // the matrix 'm':
                //        m.a_row(10)

  ~MatrixRow()    {}

                    // Assign a value to all the elements
                    // of the Matrix Row
  friend void operator  = (const MatrixRow& mr, const int val);
                      // Modify the elements in the row
  friend void operator +=  (const MatrixRow& mr, const int val);

  void operator = (const Vector & vec);    // Assign a vector to a matrix row
};

                // Construct the MatrixRow
inline MatrixRow::MatrixRow(const Matrix * mp, const int row)
  : inc(mp->nrows)
{
  matrix = mp;
  matrix->is_valid();

  register int rowind = row - matrix->row_lwb;

  if( rowind >= matrix->nrows || rowind < 0 )
    _error("Row #%d is not within the matrix",row,((*matrix).info(),0));

  MatrixRow::ptr = &(matrix->index[0][rowind]);
}

inline MatrixRow Matrix::a_row(const int rown) const
return c(this, rown)
{}


/*
 *------------------------------------------------------------------------
 *            Vector as a n*1 matrix
 */

class Vector
    : public Matrix
{
  friend class MatrixColumn;
  friend class MatrixDiag;

public:
  Vector(const short n);        // Specify a blank vector for a given
                    // dimension, with indexation at 1
  Vector(const short lwb, const short upb); // Specify a general lwb:upb vector
  Vector(const Vector& another);    // Make a vector by example

                    // Make a vector and assign init vals
  Vector(const short lwb, const short upb,  // lower and upper bounds
     double iv1, ...            // DOUBLE numbers. The last arg of
     );                 // the list must be string "END"
                    // E.g: Vector f(1,2,0.0,1.5,"END");

  ~Vector()    {}

                // Erase the old vector and create a
                // new one according to new boundaries
  void resize_to(const short n);            // Indexation @ 1
  void resize_to(const short lwb, const short upb);     // lwb:upb specif
  void resize_to(const Vector& v);            // like other vector

  REAL & operator () (const int index) const;

                // Listed below are specific vector operations
                // (unlike n*1 matrices)

                    // Status inquires
  int q_lwb(void) const        { return row_lwb; }
  int q_upb(void) const        { return nrows + row_lwb - 1; }

                    // Compute the scalar product
  friend double operator * (const Vector& v1, const Vector& v2);

                    // Vector norms
  double norm_1(void) const;               // SUM{ |v[i]| }
  double norm_2_sqr(void) const;               // SUM{ v[i]^2 }
  double norm_inf(void) const;            // MAX{ |v[i]| }

  Vector& operator = (const double val)
  { Matrix::operator =(val); return *this; }
};

                // Create a blank vector of a given size
inline Vector::Vector(const short n)
  : Matrix(n,1)
{}

                // Create a general blank vector
inline Vector::Vector(const short lwb, const short upb)
  : Matrix(lwb,upb,1,1)
{}

                // Make a vector by example
inline Vector::Vector(const Vector& another)
  : Matrix(another.q_lwb(),another.q_upb(),1,1)
{}

                // Erase the old vector body and create a
                // new one to accomodate a specified no.
                // of elements
inline void Vector::resize_to(const short n)
{
  Matrix::resize_to(n,1);
}

inline void Vector::resize_to(const short lwb, const short upb)
{
  Matrix::resize_to(lwb,upb,1,1);
}

inline void Vector::resize_to(const Vector& v)
{
  Matrix::resize_to(v.q_lwb(),v.q_upb(),1,1);
}

                // Get access to a vector element
inline REAL& Vector::operator () (const int ind) const
{
  is_valid();
  register int aind = ind - row_lwb;
  if( aind >= nelems || aind < 0 )
    _error("Requested element %d is out of Vector boundaries [%d,%d]",
       ind,row_lwb,nrows+row_lwb-1);

  return elements[aind];
}